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1. Introduction 

Quantum Chromodynamics (QCD) is a fundamental theory which describes the strong interac- 
tion of quarks and gluons — the basic blocks of matter. The knowledge of QCD at finite temperature 
and finite baryon chemical potential is essential to the understanding of a variety of phenomena. 
Exploring the QCD phase diagram, including identification of different phases and the determina- 
tion of the phase transition line is currently one of the most studied topics [jl|, g]. Guided by phe- 
nomenology and experiments, a candidate QCD phase diagram was proposed [|3j] - Fig. [T] shows 
a sketch of the conjectured phase diagram. In certain limits, in particular for large temperatures 
T or large bary on-chemical potential jx B , thermodynamics is dominated by short-distance QCD 
dynamics and the theory can be studied analytically. However, the most interesting phenomena 
lie in regions where non-perturbative features of the theory dominate. The only known systematic 
approach is the first-principle calculation via lattice QCD using Monte Carlo simulations. 




Figure 1: Conjectured QCD phase diagram 

The thermodynamics of strongly interacting matter has been studied extensively in lattice cal- 
culations at vanishing baryon chemical potential. Current lattice calculations strongly suggest that 
the transition from the low temperature hadronic phase to the high temperature phase is a continu- 
ous but rapid crossover happening in a narrow temperature interval around T c ~ 170MeV [Q, [|, 0, 
[7j> % On the other side of the phase diagram — large baryon chemical potential but very low 
temperature, a number of different model approaches [ fTl| , 12, H^j, 15, If], 17, [l8|] suggest that 



the transition in this region is strongly first order, although this argument is less robust. This first 
order phase transition is expected to become less pronounced as we lower the chemical potential 
and it should terminate at a second order phase transition point - the critical point. 

The search for the QCD critical point has attracted considerable theoretical and experimental 
attention recently. The possible existence of this point was introduced sometime ago [ ]l3[ [14]]. It 
is apparent that the location of the critical point is a key to the understanding of the QCD phase 
diagram. For experimental search of the critical point, it has been proposed to use heavy ion col- 



lisions at (RHIC) [19, 20]. The appearance of this point is closely related to hadronic fluctuations 



which may be examined by an event-by-event analysis of experimental results. The upcoming 
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RHIC lower energy scan and future LHC and FAIR (GSI) will focus on the region T > 100 MeV 
and Hb ~ — 600 MeV where the critical point is predicted to exist in theoretical models. However, 
in order to extract unambiguous signals for the QCD critical point, quantitative calculations from 
first-principle lattice QCD are indispensable. 

Remarkable progress in lattice simulations at zero baryon chemical potential has been made in 
recent years; however, simulations at non-zero chemical potential are difficult due to the complex 
nature of the fermionic determinant at non-zero chemical potential. The phase fluctuations produce 
the notorious "sign problem". The majority of current simulations are focusing on the small chem- 
ical potential region H q /T <C 1 where the "sign problem" appears to be controllable. Most of them 
are based on the grand canonical ensemble (T, }Xb as parameters). Since the existence and location 
of the critical point is still unknown, we need an algorithm which can be applied beyond the small 
chemical potential region. This is one of the motivations of our studies via the canonical ensemble 
approach [f| g g ^g]. 

We have proposed an algorithm based on the canonical partition function to alleviate the over- 



lap and fluctuation problems [E6, 27]. The method we use is computationally expensive since every 



update involves the evaluation of the fermionic determinant; however, finite baryon density simu- 
lation based on this method has been shown to be feasible [27, 28]. In this approach, we measure 
the baryon chemical potential to detect the phase transition. With the aid of the winding number 
expansion technique igg^g, a program was outlined to scan the QCD phase diagram in 
an effort to look for the critical point [33, 34]. 

In this study, we present results based on simulations using 6 3 x 4 lattices with clover fermions. 
We have run simulations using two, three and four degenerate flavors of quarks. As a benchmark, 
we first study the four flavor case. This system has a first order phase transition that extends all the 
way to zero density and the low density region has been reliable studied using staggered fermions. 
We map out the phase transition line using our method and compare it with the results of previous 
studies; we find very good agreement. For the two flavor case we present results for three different 
temperatures, the lowest one being T = 0.83r c . In this case we do not see any clear signal for 
a phase transition. While surprising, this is consistent with results from other studies [53]. The 
most interesting results we report are for the Nf = 3 case, where the evidence of a critical point 
would not only suggest the existence of such a point in the real world but also give an estimation 
of its location. We scan the phase diagram and we find signals for a first order phase transition for 
temperatures below 0.92r f . We locate the critical point at the intersection of the boundaries of the 
coexistence region. 



2. Canonical partition function 

One way to show how to build the canonical ensemble in lattice QCD is to start from the 
fugacity expansion, 

Z(V,T,li) = J £z c (V,T,k)e^ T , (2.1) 

k 

where k is the net number of quarks (number of quarks minus number of anti-quarks) and Zq is the 
canonical partition function. Using the fugacity expansion, it is easy to see that we can write the 
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canonical partition function as a Fourier transform of the grand canonical partition function, 

Z c (V,T,k) = ^J^d(l>e-^Z(V,T^)\^T- (2.2) 

after integrating over the imaginary chemical potential. 

For our study, we use the Iwasaki improved gauge action and clover fermions [36, |7|, 38]. The 



imaginary chemical potential is introduced via a U(l) phase added to the time links from the last 
time slice. 

As an illustration, we will consider the case of two degenerate flavors. After integrating out 
the fermion part, we get 

Zc(Y,T,k) = J 9\Je- s ^6et k M 2 (U), (2.3) 

where 

1 r 27t 

det k M 2 (U) = — d^e-' k UetM 2 (m,n;U)\n=^T, (2.4) 
2k Jo 

is the projected determinant with the fixed net quark number k. Using charge transformation and 
75 hermiticity property of the action, we can prove that 

Z c (V,T,k)=Z c (V,T,-k). (2.5) 

This property allows us to rewrite the partition function as 

Z c (V,T,k) = J 3HJe- s ^ U) Redet k M 2 (U). (2.6) 

The integrand is real but not necessarily positive. To compute observables based on this partition 
function we rewrite it as: 

Z c (V,T,k) = J &\Je-^ iu) det k M 2 (U) 

= J e- s ^ u) detM 2 (U)W(U)a(U), (2.7) 

where 

, . \Redet k M 2 (U)\ 

W(U) = - * \ , (2.8) 

v ; detM 2 (£/) 

and 

v ; |Redet A .M 2 ([/)| 

We separate the phase CC(U) and generate an ensemble distributed according to the positive weight 
|RedetfcM (U)\; the phase factor is reintroduced in the observable. To generate this ensemble a 
candidate configuration is "proposed" by the standard Hybrid Monte Carlo algorithm [ p9| , 50 , |4T| , 



42], then an accept/reject step is used to correct the weight. We note that the accept/reject step is 



based on the determinant ratio which alleviates the fluctuation problem [43, 23, gg] and enhances 
the acceptance rate [[2^]. 
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3. Winding number expansion 



Most of the simulation time is spent on the accept/reject step, specifically on computing deter- 
minant of the fermion matrix. On the 6 3 x 4 lattice, the matrix has 10368 x 10368 entries. Although 
it's very sparse, exact determinant calculation is numerically demanding even on this small lattice. 
Moreover, the determinant needs to be evaluated many times at every accept/reject step since we 
need to compute its Fourier transform; our original approach was to use an approximation where 
the continuous Fourier transformation is replaced by a discrete one, i.e.: 



N-l 



tet k M 2 {U) « - £ e-'^ detM(t/ ,) 2 , 



a 2Kj 



(3.1) 



It was shown that the errors introduced by this approximation are small Q34| ] for small quark num- 
bers. However, there are two problems with this approach: (i) the computation time increases 
linearly with the net quark number, (ii) although evaluating the determinant N times should theo- 
retically allow us to compute projection numbers as large as k = N/2, numerically we found that 
for large quark numbers, the Fourier components become too small to be evaluated with enough 



precision, even using double precision floating point numbers. It has been shown | J29| ] that the 
the results of the projected determinant for k larger than 20 would differ significantly for different 
choice of N, which signals a numerical instability. To see this problem, we use N = 208 to evaluate 
the fermion determinant and calculate Fourier projection using discrete Fourier transform. One 
would expect that |RedetfcM2(£/)| should decrease exponentially. We see that this is indeed the 
case for discrete Fourier transform at small quark number, as shown in Fig. || (left panel). How- 
ever, as the quark number gets close to 30, the projected determinant calculated using the discrete 
Fourier transform approximation flattens out and stops falling below the double precision limit of 
10~ 15 . This is the onset of the numerical instability. 




150 200 




Figure 2: Numerical instability of discrete Fourier transform with 208 points (Left). Comparison from 
winding number expansion method and the discrete Fourier transform with N = 16 evaluations (Right). 
Error bars are the standard deviations. 

This happens because e~' k ^> oscillates rapidly at larger k, leading to large cancellation in the 
sum over e~ lk ^> detM^^j) 2 . The accumulation of round-off errors makes it impossible to evalu- 
ate the projected determinant. This numerical challenge led us to develop a new method which 
should be free of the above numerical problem for quark numbers relevant to the phase transition 
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region in this study [p9|]. The idea of our new method is to first consider the Fourier transform 
of logdetM(£/, <p) instead of the direct Fourier transform of detM(U,(j>). Using an approxima- 
tion based on the first few components of logdetM(£/, 0), we can then analytically compute the 
projected determinant. The efficacy of the method can be traced to the fact that the Fourier compo- 
nents of logdetM(£/, 0) are the number of terms in the expansion which characterizes the number 
of quark loops wrapping around the time boundary. They are exponentially smaller with increasing 
winding numbers. This is why we can approximate the exponent of the determinant very accurately 
with a few terms which, in turn, allows us to evaluate the Fourier components of the determinant 
precisely. 

To see how this works, we look at the hopping expansion of the logdetM(£/, <f>). We start by 
writing the determinant in terms of the trace log of the quark matrix 

detM(t/,0) =exp(TrlogM(t/,0)). (3.2) 

It is well known that Tr logM corresponds to a sum of quark loops. We separate all these loops 
in classes depending on how many times they wrap around the lattice in the temporal direction. We 
have then 

TrlogM(l/,4>) = £ L(U,$) (3.3) 

loops 

= A (U) + [^e^W n {U)+e- in ^(U)l 

n 

where n is the winding number of quark loops wrapping around the time direction and W n is the 
weight associated with all these loops with winding number n. Ao(U) is the contribution with zero 
winding number. Eq. ( [0|) can be re-written as 

TrlogM(t/,0) =A (U) + [^e^W n (U) + e- in ^(U)] (3.4) 

n 

= A (U) + ^A n cos(n<j) + S n ), 



where A„ = 2\W„ | and <5„ = arg(W„) are independent of 0. Using Eq. ( p?4 ) and Eq. (fTj) we get 

detM(U,(j>) = e Ao+A l cos(t+8 l )+A 2 cos(2$+S 2 )+ _ (3 5) 

The Fourier transform of the first order in the expansion can now be computed analytically; 
we have 

r n ftyg-a^+A, = e Mh h i Ax \ (3 6) 

Jo 2% 

where 4 is the Bessel function of the first kind with rank k. 

For higher orders in the winding number expansion, we compute the Fourier transform based 
on the truncated Taylor expansion: 

2% 

f-2n A ,k 00 <*> A n k 

' a( P —ikd> „An+Ai cosf(4+5i ) TT V A k 







[■lit A A °° 00 A K 

/ ^ e -^/o+A 1 co S (0+5 1 ) -Q £ lk. cos (jk0 + 
27C k=? „,.=n n k- 



k=2n k =0 ■ 

cool* (A i) +c+oiIfc+i(Ai) + c_oiI*-i(Ai) + c+02l/t+2(Ai) + ... (3.7) 
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The projected determinant is written in terms of the linear combination of Bessel functions, the 



coefficients c can be easily computed analytically. Using Eq. ( |3.7[ ) and the recursion relation for 
the Bessel function, (A) = ^I&(A) +Ijt+i (A), the winding number expansion method (WNEM) 
can be extended to higher orders. 

As we mentioned before, the efficacy of the method rests on the fact that we can get a very 
good approximation of the exponent using only a few terms in the Fourier expansion. However, 



the evaluation in Eq. (3.7) is based on a mixture of analytical integration and a truncated Taylor 



expansion. Since the error of the truncated Taylor expansion is not well controlled, to produce 



our final results, we decided to employ a method which can evaluate Eq. ( p.l\ ) more precisely and 
thus correct the necessary errors introduced by the Taylor approximation. The methodology can be 
viewed as a form of reweighting and the approximation level is then tuned to produce errors that 
are smaller than the statistical errors of the final result. The strategy is the following: 

• During the procedure of generating ensembles, we compute the determinant for 16 phases 
and we keep only 6 winding loops in the TrlogM expansion which has been shown to be 



precisely enough for the approximation of TrlogM [|29|]. The evaluation of the projected 
determinant is done by truncating the second and third winding loops terms to their tenth 
order, and first order for the remaining terms. We denote the approximated determinant 
using truncated Taylor series in Eq. (|3.7|) by det^M. 



Once the ensembles are generated, we evaluate the determinant of each configuration for 
No phases which admits projection up to support as large us Nq/2 winding loops reliably. 



Instead of evaluating the integral Eq (3.7) by the truncated Taylor series, we evaluate it 
numerically through a multi-precision library (GMP library) by setting precision as 512 digits 
after decimal point. In order to choose No, we include the higher order of winding loops until 
the projected determinant converges, the minimal number of winding loops included will be 
equal to No/ 2. 

Since the ensemble generated in the Monte Carlo process involves an approximation of the 
projected determinant and is not the target ensemble, we correct for this in the observable 
with the following reweighting 

(0) = ^ (3.8) 
{a r )o 

where ()o denotes the average over the ensemble generated with the approximated determi- 
nant and 

c,= £^ (3.9) 
det^M 

is the reweighting phase where detkM is from the exact Fourier transform without Taylor 
expansion. If the approximation due to Taylor expansion is very good then (a r )o will be 
very close to 1; if the difference between the truncated determinant and the exact one is 
much smaller than the error bars on (O) then the reweighting step is not really needed and 
we can consider the approximation exact. On the other hand, if the average is significantly 
different from 1 we need to do reweighting. The worst case is when the approximation is so 
bad as to introduce a sign problem ((a,-)o overlaps with zero within the error bars). In that 
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case the approximation is invalid since the generated ensemble cannot be used to compute 
any observable. 

The detailed discussion of this methodology will be presented somewhere By employing 

the above strategy, even if the approximation is not perfect, it can be corrected by a reweighting 
after the ensemble is generated. The only possible issue is to make sure that the reweighting phase 
doesn't introduce too much noise in the observable; if that is the case, the approximation needs to 
be improved and a new ensemble needs to be generated. 

Before we conclude this section, we would like to compare the merits of WNEM with the those 
of discrete Fourier transform. We compute the values of the projected determinant determined 
using the discrete Fourier transform with only N = 16 in Eq. ( p.l[ ). In Fig. || (right panel), these 
results are labeled discrete (16). We see that this approximation is only valid up to k = 6. This is 
to be compared to WNEM (16) which takes the same computational time and yet does not suffer 
from this problem and, as expected, the evaluated determinant continues to decrease as the quark 
number is increased. It is this projected determinant evaluated from WNEM that allows us to scan 
a wide range of densities on the QCD phase diagram. 



4. Results 

4.1 Baryon chemical potential 

Before starting the discussion on the QCD phase diagram, we will first present the "tool" we 
used to determine the phase transition in the canonical ensemble: the baryon chemical potential. 
We measure the chemical potential and plot it as a function of the net quark number k. Due to 
the contribution of the surface tension to the free energy at finite volume, the chemical potential 
in the mixed phase region has an "S-shape" structure. We scan the phase diagram looking for this 
S-shape a structure to signal a first order phase transition. 

In the canonical ensemble, the baryon chemical potential is measured by taking the difference 
of the free energy after adding one baryon, i.e. 



F(n B + \)-F(n B ) = 1 Zc(3n B + 3) = 1 (/(£/))„ 
(n B + l)-n B P n Z c (3n B ) j3 {a(U)) 



where 



, , Redet 3nR M n f(U) 

a(U) = i — h4r, and (4.2) 

V j |Redet 3 , Js M"/ (£/)[' 

Rede^M^l/) 

7{U) |Redet 3 „ B M"/([/)| ( V 

are phase factors measured in the ensemble with n B baryon number where () c in Eq. (4.1) stands 
for the average over the ensemble generated with the measure |Redet3„ a M"^ (U)\. 

4.2 Phase diagram for N f = 4 and N f = 2 

We turn now our attention to the QCD phase diagram. Although it has been speculated for 
quite a some time that the QCD phase diagram may look like Fig [I], unfortunately, there is little 
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solid evidence to to support this scenario. The finite temperature studies at zero baryon density 
have shown that the transition from hadronic matter to quark gluon plasma is a rapid cross-over. 
Lattice QCD evidence regarding the existence and location of the second order transition point will 
help tremendously in sharpening our understanding of the QCD phase diagram. 

Before we present our results, we would like to point out the difference between the phase 
diagram in the grand canonical ensemble and the one in the canonical ensemble. We plot the 
canonical partition function phase diagram in Fig. |[ In terms of T and p, the first order phase 
transition line becomes a phase coexistence region which has two boundaries that separates it from 
the pure phases. As we approach the critical point, the two boundaries will eventually meet at that 
point. We determine the boundaries of the coexistence phase at different temperatures and locate 
the critical point by extrapolation. 

As a first attempt, we carry out simulations with four and two degenerated quark flavors. The 
reason we start with two and four flavors is the following: first of all, it is easier to simulate an 
even number of flavors with HMC [40]. Secondly, in the four-flavor case, the first order phase 
transition line extends all the way to zero density. This provides a testing ground for the algorithm 
since there are reliable results mapping the transition line in the small chemical potential region. 
Thirdly, the two-flavor case is expected to have a phase diagram similar to that of real QCD. The 
expected phase diagrams for two and four flavor cases are shown in Fig. [|. For this reason, we 
will use the four flavor simulation as a benchmark to demonstrate the methodology of determining 
the phase boundaries and locating the critical point before tackling the more realistic case. The 
expected phase diagrams for two and four flavor cases are shown in Fig. |3[ 

Given the first order transition at zero chemical potential for Nf = 4, it is natural to ask whether 
this first order is preserved when we switch on the chemical potential. As we mentioned before, 
the first order phase transition is reflected in an "S-shape" structure in the vs. p plot (or Hb vs. 
rig plot as we will present). To proceed in this way, we decide to scan the phase diagram by fixing 
the temperature below T c while varying n%. 

We scan the phase diagram by fixing the temperature below T c while varying wg. Once we 
cross into the coexistence region, in finite volume, the non-zero contribution from the surface 
tension causes the appearance of a "double-well" effective potential [^] which leads to an S- 
shaped behavior in the chemical potential plot. In the thermodynamic limit, the surface tension 
contribution goes away since it is a surface term and the free energy scales with the volume; the 



9 



Study of QCD critical point using the canonical ensemble method 



Anyi Li 



chemical potential will then stay constant in the mixed-phase region. 



T=0.89T r 




* * 
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Figure 4: Nf — 4 phase scan along with Maxwell constructions 

We present our results for Nf = 4 in Fig. |] for three different temperature below T c . To identify 
the boundaries of the coexistence region and the critical value for the baryon chemical potential, we 



rely on the "Maxwell construction" [25]: we select several points in the S-shape region and fit them 
with a third-order polynomial. We determine the critical chemical potential pL c and the boundary 
points pi and P2 by requiring the two areas between the fitted curve and a constant horizontal line 
which intersects the curve at pi and P2 to be equal. For all reasonable fits, we find the values of 
the boundary points, p\ and P2, and the value of the critical chemical potential (/i c ) to be fairly 
insensitive to our choice of the fit function or fit region. The simple third order polynomial fit is 
sufficient. 

We perform the Maxwell constructions for the three temperatures we studied: 0.89, 0.91 and 
0.95T C ; the results are presented in Fig. |J Having determined the pi and P2 for three different 
temperatures, we plot boundaries of the coexistence region by a simple extrapolation. Although 
there is no "critical point" for the Nf = 4 case, it is expected that the two coexistence phase bound- 
aries should cross at zero chemical potential and T = T c . To determine the crossing point, we fit 
the boundary lines using a simple even polynomial in baryon density 

^ = 1 - a(N f ,m q )(Vp) 2 + O ((V p) 4 ) (4.4) 

to do the extrapolation. The phase boundaries and their extrapolation are plotted in Fig. |5[ We find 
the intersection point at T c (p)/T c (0) = 1.01(5) and p = 0.05(10)fm~ 3 which is consistent with 
our expectation: pL = and T = T c . 

We would like to point out that the critical temperature was determined using a different set of 
simulations. We run simulations at zero baryon density and we monitored the Polyakov loop sus- 
ceptibility as we increased the temperature. These simulations were performed using the standard 
HMC algorithm and we used two lattice sizes: 6 3 x 4 and 10 3 x 4. The critical temperature was 
determined using the peak of the Polyakov loop susceptibility and we also checked that the suscep- 
tibility peak scales with the volume as expected. The fact that the intersection point is close to the 
critical temperature determined from a set of simulations using a different methodology constitutes 
a convincing cross-check of the correctness of our method. 

We compare our results to the study of analytic continuation from the imaginary chemical 



potential with stagger fermions [45]. The results are presented in Fig. n (left panel) - we see that 



the agreement is quite good in spite of the fact that the quark mass used in the staggered fermion 
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Figure 5: Phase boundaries in the temperature vs. density plot for Nf = 4. 



study, m n 300 — 400 MeV, is significantly different from the one we use. This is not surprising 
since the phase transition curve for Nf = 4 is fairly independent of the quark mass p6[|. As a check 
that the reweighting does not change drastically the observable results, we also plot the results with 
no reweighting which should contain systematic errors from Taylor expansion. Nevertheless, we 
see that the errors of the chemical potential from non-reweighted results overlap with those from 
the reweighted ones and the extrapolated crossing point of the two boundaries also agree with T c at 
/i = albeit with larger errors. We also compare our results to ones from a canonical ensemble 
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Figure 6: Phase transition line for Nf — 4 in the T, jj. plane. 

study that uses staggered fermions [^] (Fig. ^| (right panel)). We find that the results are consistent 
with each other; however, our coexistence region is narrower. The discrepancy could come from 
the difference in the quark mass (m % « 800 MeV for our study compared to m n th 300 MeV for the 
staggered study), and/or the fact that we use a different fermion formulation. 

As a sanity check, we examine the seriousness of the potential sign problem. The average 
sign in Eq.(3.S) is plotted in Fig. [7| (left panel). We see that they range mostly from 0.6 to 0.2. 
Except for the point at T = 0.89T C and jib = 9 which has a 1.5 sigma away from zero, the others 
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Figure 7: Left: Average sign of Nf — 4 in terms of ng- Right: Average sign of Nf — 2 in terms of «g 

all have more than 3 sigmas above zero. In view of the fact that our results agree with those from 
the imaginary chemical potential study and the extrapolated boundaries meet at the expected T c at 
jj. = 0, we believe that the sign fluctuation in the reweighting is not a problem. 
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Figure 8: Baryon chemical potential for Nf = 2. 

For the Nf = 2 case, the phase diagram (see Fig. ^ is expected to be similar to the conjectured 
real QCD phase diagram. For nonzero quark mass, it also shows a crossover behavior at vanishing 
chemical potential, a critical point as well as a first order phase transition line. For this system, we 
scan three temperatures below T c : 0.83, 0.87 and 0.92T C . The results are presented in Fig. ||: Given 
our limited statistics of 2000 configurations and the fact that the sign problem starts to set in below 
0.83T C , we do not observe a clear signal for the S-shape structure within the scanned temperature 
and density range. This may be due to the fact that the transition is milder than the sensitivity 
of our data, or it may be due to the fact that the transition occurs at lower temperatures than the 
ones we studied - there is at least one study that claims that the critical point for Nf = 2 occurs at 



temperatures below 0.8T C []35|]. However, in view of the sign problem we encounter below 0.83T C , 
the claim needs to be scrutinized. 

We show the average sign in Eq. (3^) in Fig. ^ (right panel). Again, we see that except for the 



case at the higher «g and lower temperature, the average signs for most of the cases are more than 
three sigmas above zero. 

4.3 Phase diagram for Nf = 3 

As we see from the above studies, the canonical approach at finite density is feasible and has 
been checked in the Nf = 4 case. We would like to extend this study to the real world which 
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contains two light quarks and one strange quark. However, simulating light quark masses requires 
larger volumes than the ones used in the present study. As a first attempt, we shall use the quark 
mass around physical strange quark mass as our input and apply the canonical ensemble method 
to the three degenerated flavor case which is expected to have similar features to the real QCD 
case. We are looking for the existence of a first order phase transition and its critical point at this 
stage, which may serve as an indication for the existence and approximate location of such a point 
in the real world. Following the methodology for two and four flavors, we fix the temperature at 
four different values and scan for the phase diagram by varying n^. The resulting baryon chemical 
potential is plotted in Fig. |9[ 
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Figure 9: Nf = 3 phase scanning along with Maxwell constructions 

The presence of the S-shape structure at three temperatures below 0.92T C is a signal for a first 
order phase transition. Determination of the coexistence phase boundaries is carried out via the 
Maxwell construction. We plot the results from the Maxwell constructions in Fig. |9| 

In Fig. |l^ (left panel) we plot the boundary points as determined by the Maxwell construction. 
The phase boundaries intersect at the critical point which is located at T/T c = 0.94(3) and pe = 
1.82(7)fm~ 3 . In the right panel of Fig. [E], we plot the critical chemical potential as a function of 
temperature. On this graph the critical point is located at Te = 0.94(3)T C and \Le = 3.01(12)T C . 

Now that we have determined the phase boundaries in the temperature-density plot, we could 
map out the phase boundary in the temperature-chemical potential plot. The later is the usual phase 



diagram studied in the grand canonical ensemble. It is shown in Fig. 10 (right panel). 

From our current simulations, we have found a first order phase transition and determined the 
location of the critical point for the three degenerate flavor case at intermediate quark mass. At the 
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Figure 10: Phase boundaries in the temperature vs. density plot for Nf = 3 (left). Transition line in the 
temperature vs. chemical potential plot (right) 
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Figure 11: Average sign of the Nf = 3 case for different ng and T (left: with reweighting; right: without 
re weighing). 



critical point, the critical temperature is determined as Te = 0.94(3)T C and the baryon chemical 
potential is He = 3.01(12)T C . Given Te = 0.94(3)T C , we are supposed to see an S-shaped behavior 
at T = 0.92T C . However, we cannot observe a clear signal of it, it may be due to the fact that this 
temperature is too close to Te- 

We show the average sign from reweighting (Eq. (3.S)) in Fig. |ll] (left panel). Again, we see 
that most of them are more than 3 sigmas above zero. The only exceptions are those points at 
hb > 14 which are beyond the coexistence phase region. We also plot the average sign with no 



reweighting (Eq. (|2.9|)) in Fig. 1 1 (right panel). As expected, we find that they are in general larger 
in value than their respective ones with reweighting. 

Before we conclude this section, we would like to discuss a potential contradiction between 
our results and those of de Forcrand and Philipsen |48| , fi^ ]. These authors used 2+1 and 3 
flavors staggered fermions and a Taylor expansion in H q /T to study the curvature of the critical 
surface at very light quark masses close to \i q = surface. After expanding to {pL q /T) A , they find 
that the critical surface bends such that the first order region shrinks at higher quark masses. This 
suggests an exotic scenario where there is no critical point at finite chemical potential [ 47 , 48 , [4^] . 
This conclusion seems to contradict our results. However, since this conclusion is derived from 
studies at very low chemical potential, it is possible that the critical surface bends back at larger 
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chemical potentials and the critical point reappears. In fact, an NJL model study seems to indicate 



that this scenario is actually realized [ |50| ] 



5. Summary 

In this study, we show that the canonical ensemble can be used to investigate the phase diagram 
at finite temperature and nonzero baryon density, at least on small lattices. In order to scan the QCD 
phase diagram at relatively large quark number (e.g. k > 20), we developed the winding number 
expansion method to calculate the the projected determinant. It greatly expands our ability of 
scanning a broad region of the phase space and allows us to reach the coexistence phase. 

We presented our results from simulations using two and four degenerate flavors of quarks. 
At finite volume, the first order phase transition at finite density appears as an S-shaped structure 
in the plot of chemical potential versus baryon number (baryon density). To determine the phase 
boundaries we use the Maxwell construction. For Nf = 4, the first order phase transition is observed 
at three different temperatures below the transition temperature at zero chemical potential. The 
phase boundary separating the hadron gas phase from the coexistence phase and the boundary 
between the plasma phase and the coexistence phase eventually cross at one point. For the four 
flavor case, this point should be the first order transition point at zero chemical potential. We have 
checked this numerically and found that the intersection of the boundaries indeed coincides with 
the the critical temperature at /I = within errors. Furthermore, our results are consistent with 
those of the imaginary chemical potential approach in the grand canonical ensemble [B^]. This is 
a powerful check for our method. For the two flavor case, we do not see any clear signal of a first 
order phase transition for temperature as low as 0.83T C . 

The most interesting study is the three flavor case. We found signals of a first order phase tran- 
sitions at three different temperatures and determined the phase boundaries for each temperature. 
The critical point located at T = 0.94(3)T C and jU# = 3.01(12)T C is obtained by an extrapolation. 

We should point out that the present simulations are conducted on a small volumes, and with 
relatively large quark masses and a coarse lattice spacing (a ~ 0.3 fm). We plan to run simulations 
at lower quark masses using nHYP clover fermions. 
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